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Abstract. The possibility that some Markarian objects (e.g. Mkn 501, Mkn 421 and Mkn 766) host massive 
binary black hole systems with eccentric orbits at their centers has been considered. These systems could be 
sources of gravitational radiation for space-based gravitational wave interferometers like LISA and ASTROD. In 
the framework of the Lincoln - Will approximation we simulate coalescence of such systems, calculate gravitational 
wave templates and discuss parameters of these binary black hole systems corresponding to the facilities of LISA 
and ASTROD. We discuss also the possibility to extract information about parameters of the binary black hole 
systems (masses, of components, distances between them, eccentricity and orbit inclination angle with respect to 
line of sight) from future gravitational wave measurements. 
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1. Introduction 

BL Lacertae objects, also known as Markarian objects 
(hereafter Mkns), belong to the class of active galaxies 
1 according to the well-esta blished unified model on radio- 
loud active galactic nuclei l|Urrv fc Padovanill995|) . These 
objects are thought to be dominated by relativistic jets 
seen at small angles to the line of sight. 

Until now, several astrophysical phenomena have 
been attributed to binary black holes, like precession 
of jet s (BcEclman, Blandfo rd fc ReesHl980j) . misalignment 
llConwav fc Wroblel I199.3T periodic outburst activity in 
the Quasar OJ 287 llSijlanpaa. Haarala. Valtonen et alJ 
Il988l iLehto fc ValtonerT TT996) and pre cession of t he ac- 
cretion disk under gravitational torque llKatzlll997lh 

It has been recently observed that some Mkns show a 
periodic behavior in the radio, optical, X-ray and 7-ray 
light curves that is possibly related to the presence of a 
massive binary black hole with a jet along the line of sight 
or interacting with an accretion disk ( Yu 2002). Therefore, 
the search for light curve variability, mainly in X-ray and 
7-ray wavelengths, can be considered as a method to probe 
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the existence of a massive binary black hole in the center 
of a galaxy. 

A question naturally arises: how can a binary system of 
massive black holes be formed? The answer to this prob- 
lem can be found in the framewo rk of the hierarchical vi- 
sion of the universe l|Whitelll997jh for example if massive 
black hole systems form as a result of merging processes 
between galaxies, eac h of them may contain in the center 
a massive black hole jReeslll984l Kormendv fc Richstond 
Il995l iRichstone. Aihar fc Bender! 1 1998|) . Recent observa- 
tional sign atures 1 of such hypothesis were discussed, for 
example bv lYu fc LiJ l|200lj) who analyzed typical features 
of Fe K a line shapes. 

At least three Mkns (i.e. Mkn 501, Mkn 421 and Mkn 
766) are particularly well studied at high energies, reveal- 
ing a possible periodic behavior in their light curves. 

Mkn 501, at z — 0.034, shows a clear well-correlated 23 
day periodicity in X-ray and TeV energy bands with an 
observed TeV flux ratio / — 8 between the maximum and 
minimum of the signal l|Protheroe. Bhat fc Fleurvl 

1 The gravitational wave spectrum from coalescing massive 
black hole binaries formed b y merging processes o f their host 
galaxies has been studied bv lJaffe fc Backer! <|2003T) . 
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IHavashida. Hirasawa fc Ishikawal 

19981 IKranich. DeJager. Kestel et all " | l999f 

Nishikawa. Havashi. Chamoto et alJ Il999l: 
Kranich. DeJager. Kestel et al while evi- 

dence for cor relations in the optical U-band is 
rather weak ( [Catane^e^^radbarr^Beshr^etal.1 Il997t 
lAharonian. Akhperianian. Barrio et al.l ll999altJ ). It has 
been also suggested that the complex morphology of 
the jet and the peculiar behavior of its spectral energy 
distribution are probably rel ated to the presence of 
a massive binary black hole llConwav fc Wroblel 119951 
IVillata fc Raiteril ll999") at the center of Mkn 501. 

Mkn 421, at z = 0.031, is the brightest BL Lacertae 
object at X-ray and UV wavelengths and it is the 
first extragalactic source dis co vered at TeV energies 
l|Punch. Akerlof. Cawlev et all Il992f) . This nearby 
source, whi ch has been recently observed by th e XMM- 
Newton ijBrinkmann. Sembav. Griffiths et alJ l200ll) 
and by Beppo-SAX raschi. Fossati. Tavecchio et alJ 
Il999l) satellites, shows remarkable X-ray variabil- 
ity correlated with strong activit y at TeV energies 
ijGeorge. Warwick fe Bromagj Il988l) on a time-sca le of 
~ 10 4 s ijMaraschi. Fossati. Tavecchio et al.1 Il999l) and 
with a flux ratio / — 2. 

X-ray observations of the nearby Mkn 766, at z = 
0.01 3, have been performed by the XM M-Newton satel- 
lite l|Boller. Keil. Trumper et aTl l200lj) . These observa- 
tions have revealed the presence of a strong X-ray peri- 
odic signal with frequency ~ 2.4 x 10 -4 Hz and flux ratio 
/-1.3. 

Based on the assumption that the periodic behavior of 
the observed light curve for Mkn 501 is related to the pres- 
ence of a binary system of black holes (one of which emits 
a jet moving with Lorentz factor 7b) on circular orbits, 
iRieger fc Mannheim I l|2000ll have proposed a method to 
determine the physical parameters of the binary system 
from the observed quantities, i. e. signal periodicity, flux 
ratio between maximum and minimum signal and power 
law spectral index. However, binary black holes might 
be on eccentric orbits and eccentric ity values up t o 0.8- 
0.9 are not necessarily too extreme ljFitchetdll987h . This 
could happen frequently if binary systems of black holes at 
galaxy centers form generally as a consequence of merging 
processes between their host galaxies. 

Of course, due to gravitational wave emission, orbits 
tend, as a first approximation, to circularize but this hap- 
pens within a time-scale of the same order of magnitud e 
as the merging time-scale freterd ll 964 iFitchettl ll 987K 
Therefore, if a massive binary black hole is found at the 
center of a galaxy, it may happen that the constituting 
black holes are still on eccentric orbits, in which case the 
method proposed bv lRieger fc Mannheim I l|2000h does not 
hold. 

Hence, iDe Paolis. Ingrosso fc Nutital l|2002j) studied 
the massive black hole binary system possible in the center 
of some Mkn objects by assuming more general elliptical 
orbits and considered the orbit eccentricity e, the binary 
separation a and Lorentz factor 7b as free fit parameters 



used to determine the masses (Mi and M2) of the two 
black holes from the observed X-ray periodicity. Once the 
orbital parameters of the MBH (i.e. Massive Black Hole) 
binaries are known, the values of the obtained orbital sep- 
aration, eccentricity and MBH masses are considered as 
initial conditions in the time evolution of the binary sys- 
tems. 

In the present paper we study the evolution of the sys- 
tem and to determine the gravitational wave waveforms, 
i.e. the amplitude of the metric perturbation as a function 
of time. In doing this, we simul ate the evolution of bi nary 
black hole systems by using the lLincoln fc Willi il99(t) ap- 
proximation and calculate gravitational wave templates 
without any assumptions about the evolution of our sys- 
tem on quasi-circular orbits. We also study the detectabil- 
ity of the emitted gravitational waves by the n ext gener- 
ation of space-based i nterferometers like L ISA ( Rcinhard 
l200(t and ASTROD llWu. Xu fc Ni1l200(i E^ l2002l> . 

The paper is structured as follows: in Section 2 we 
show how to determine the massive black hole binary sys- 
tem parameters starting from the observed X-ray peri- 
odicity toward the considered Mkn objects. In Section 3, 
the model we use to simulate the binary system evolution 
in the post 5/,2 -Newtonian approximation is reviewed. In 
Section 4 we describe our results about the evolution of 
binary system, profile gravitational wave templates, typ- 
ical times for the evolution of our binary systems before 
fitting into the LISA frequency band. Finally, in Section 
5 we draw some conclusions. 



2. Initial conditions for a coalescing binary 
massive black hole system 

As discussed above, there exist a number of Mkn objects 
(like Mkn 501, Mkn 421 and Mkn 766) that show periodic 
activity in the radio, optical, X-ray and 7-ray light curves. 
It has been recently proposed that the observed periodic- 
ity is possibly related to the presence of a massive binary 
black hole creating a jet along the line of sight. 

Following Rieger fc Mannheim I l|2000fl . we assume 
that the observed signal periodicity has a geometrical ori- 
gin, being a consequence of a Doppler-shifted modulation. 
It is therefore possible to relate the observed signal period 
P bs to the Keplerian orbital period Pk by 



Pobs = (1 + z) ( 1 - — cosi ) Pk , 



(1) 



where i ~ l/7b and v z is the typical jet v elocity assumed 
to be u z ~ c(l - l/7b) 1/2 l|Spadal 11999]) . Here, we rely 
on the assumption that the lighter black hole in the bi- 
nary system is emitting a jet which is moving toward the 
observer with Lorentz factor 7b. 

The observed flux modulation due to Doppler boosting 
can be written as 



S(u) = <5 3+Q S» 



(2) 
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Fig. 1. The secondary mass to (in the text Mi) is shown as 
a function of the primary mass M (in the text M2) for the 
black hole binaries at the center of Mkn 501, Mkn 421 and 
Mkn 766 (from the upper to the bottom panel). The thick 
and thin lines represent the conditions expressed in eqs. 
(J7J with Lorentz factor 7b = 10 and 7b = 15, respectively. 
The binary semi-major axes a are set to 5 x 10 16 cm, 
2 x 10 14 cm and 1 x 10 13 cm for Mkn 501, Mkn 421 and 
Mkn 766. The orbital eccentricity in all cases is set equal 
to 0.5. The intersection between lines corresponding to the 
same Lorentz factor gives the masses of the black holes in 
the considered binary system. 



where a is the source spectral index 2 and the Doppler 
factor 8 is given by 



6 = 



1 — (v z cos i + vis sin i) /c 



(3) 



Here, v z is the jet velocity, i is the inclination angle be- 
tween the jet axis and the line of sight and v\ s the compo- 
nent of the less massive black hole velocity along the line 
of sight. 



Values of the power law index a are found to be 
1.2, 1.7 and 2.11 for Mkn 501, Mkn 421 and Mkn 766 
respec tive ly. For more details see iRieeer fc Man nheim I 
feOOOf l . iGuainazzi. Vacanti. Malizia et al.l (TT99sf) and 



iBoller. Keil. Triimper et al.l (l200lD . 




Fig. 2. Primary and secondary black hole masses are plot- 
ted as a function of the orbit eccentricity e for the binary 
system in Mkn 501. Dotted lines correspond to Lorentz 
factor 7b = 10 while solid lines correspond to 7b = 30. 
The semi-major axes is a = 8 x 10 16 cm. Similar plots can 
be obtained for the Mkn 421 and Mkn 766. 



Depending on the position of Mi along its orbit, the 
velocity v\ s ranges between a minimum and a maximum 
value corresponding, through eq. 10), to the two extremal 
values for the Doppler factor given by 



1 _ 



(v 2 



COS I 



27r_R a sin i 
Pk(l-e 2 )^ 2 



(4) 



K 2 +Q 



27r_R a sin i \ r —l 
Pk(l-e 2 ) 1 / 2 1 ' 



where _R a = M2a/(Mi + M2), being a and e the binary 
semi-major axis and the orbit eccentricity, respectively. 
From eq. , the observed maximum and minimum fluxes 
modulated by the Doppler effect turn out to be 



(5) 



so that one obtains the condition 5 ma x/<5min — f x ^ 3+a \ 
where / = S max (v) / S min (v) is the observed maximum to 
minimum flux ratio. 

By using eqs. © and © we have 



cPk f 1/3+a -1/1 



2vr /V3+a + 1 ^sini c 



coU I (1 - e 2 )^ . (6) 



Finally, eqs. Q and © can be rewrit- 
ten in the following form (for more details see 
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- 1 1 COS 2 
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(7) 



Mi + M 2 = 



27r(l + z) \1 cos i 



a 
G 



In these equations 7b, e and a have to be considered as 
free model parameters to be determined by the observed 
data (P bs, / and a) of the three Mkns of interest. In this 
way, we are able to estimate the masses of the black holes 
supposed to be composing each of the MBH binaries. For 
the three considered Mkn objects, the secondary mass Mi 
as a function of the primary one M 2 is shown in Fig. ^ 
The intersection between lines corresponding to the same 
Lorentz factor (7b = 10 for the thick lines and 7b = 15 for 
the thin lines) gives the masses of the black holes in the 
considered binary system. 

In Fig. |3 (for Mkn 501) the primary and secondary 
black hole masses as a function of the orbit eccentricity 
(7b = 10, 30 and a = 8x 10 16 cm) are shown. In this case, 
the black hole masses are in the range 10 6 — 10 9 M Q . 

The previous method allows us to determine the or- 
bital parameters of the massive binary black hole possibly 
in the center of some Mkn objects. The obtained values for 
the semi-major axes a and eccentricity e obviously changes 
in time as a consequence of emission of gravitation al radi- 
ation. However, as shown bv lLincoln fc Willi fl 99m. when 
studying binary systems non-negligible relativistic correc- 
tions in the equations of motion have to be considered. 
It follows that in the simulation of the evolution of the 
binary systems, the obtained values for the semi-major 
axes and eccentricity have to be considered as the initial 
conditions of the problem, i.e. ai and ei, respectively. 

3. Model for a coalescing binary system evolution 

To analyze the binary black hole evolution we use 
a (post) 5//2 -Newtoni an approximation developed by 

ILincoln fc Willi lll990l). 3 

Earlier. IGrishchuk fc Koneikinl l|l98(i() derived similar 
equations of motion for isolated bodies in a (post) 5 / 2 - 
Newtonian approximation taking into account the effects 
of gravitational radiation friction which lead to orb ital 
shrink . In studying these effects iGrishchuk fc Koneikinl 
( 198|J) derived their equations using oscula ting elements 
like those reported i n | Lincoln fc W ill (1990). Here we use 
the ILincoln fc Willi jl99(t approach since it is very conve- 
nient for our computational aims. Of course, there are a 



3 One can see also IZakharovl dl996l. Il99l Il998l) where this 
approximation was used to analyze the evolution of a binary 
neutron star system in which the lower mass object is losing 
matter since it fills its Roche lobe. 



lot of other more precise approaches which may allow one 
to describe the final stages of evo lution of binary black 



hole sys tems (se e, for example, | Buonanno fc Damour 
1 19991 120001); iFTziev fc Todorovl <l200ll); iDamour 



200 if: iBuonannol <l2002l): iMora fc Willi fcOO 



[ Gourhoulhon. Grandclement fe B onazzoT al ll2002alli 
lO'Shaughnessvl l)2003|) : iBuonanno. Chen fc Valisnerl 
l|2003]i and references therein). Even the transition 
from inspir al stage to plunge for circular orbits was 
disc ussed bvlOri fc Thornel I2OOO) and for eccentric orbits 



bv lO'Shaughnessvl (120 03V Idampcdaki s fc Kennefickl 
ll2002|) ; iGlampedakis. Hughes fc Kennefickl l)2002|) used 
a Kerr metric approximation to analyze the evolution 
of a binary black hole with a small mass ratio, namely 
m/M ~ 10~ 4 — 10 -6 , where m and M are the masses 
of the captured body and the central black hole, re- 
spectively. Rapidly spinning massive black hole binary 
systems as po ssible sources for LISA were considered by 
IVecchiol l)2003h in the framework of the post 15 -Newtonian 
approximation and it was found that the black hole spin 
could be a very essential factor and could drastically 
change the signature of the gravitational wave signal. 
However, in this paper we decided to use the Lincoln - 
Will approximation since, as it was mentioned before, 5/2 
is the first order of the post-Newtonian approximation 
where one could introduce the gravitational radiation 
friction in a self-consistent way. 

Moreover, usually these binary black hole systems are 
so far from the plunge stage of their evolution that it 
can be simulated by using the lowest self-consistent post- 
Newtonian approximation. Thus, one could use more sim- 
ple tr e atments like t hat proposed bv IPeters fc Mathewsl 
lll963|) : IPetCTsl l)l964) . However, in the framework of these 
approximations, the motion and the time evolution of 
the binary system cannot be studied in detail since some 
relativistic e ffects, such as the "perihelion shift" (see, 
for ex ample, lEhlers et alJ l)l976(k IGrishchuk fc Kopeikinl 
; lOMii for d e tails), are not accounted for. In addition, 
Ehl ers et al.l l|l976|) have some doubts about the accu- 
racy of such methods s ince approaches sim ilar to the 
IPeters fc Mathewsl (|l963j) and IPeters! l)l964|) approxima- 
tions are mathematically inconsistent and may lead to 
errors of the same order as the effects being considered. 
Thus, for a more accurate treat ment of the problem we 
follow the ILincoln fc Willi (|l990^> approach, to which we 
refer for more details. 

The orbital motion of two massive black holes mov- 
ing around the common center of mass is strongly influ- 
enced by gravitational radiation losses. Viewing this as a 
Newtonian orbital motion with perturbation suggests use 
of the osculating orbital element method taken from ce- 
lestial mechanics. 

The basic scenario is the following: at any time one can 
find a Keplerian orbit which is tangent to the true orbit. 
This means that both the position and the velocity of the 
particle on the true orbit coincide with the position and 
the velocity of the tangent Keplerian orbit at the consid- 
ered time. Of course, at a subsequent instant the actual 
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a) b) 

Fig. 3. For Mkn 421, the evolution of a MBH binary system is shown. The orbital parameters reported in each panel 
allow to reproduce the X-ray light curve periodicity. The evolutions of both the MBH binary system eccentricity a) 
and separation b) are given as a function of the number of revolutions $/27r. The initial parameter pi = a;(l — e?) is 
given in unit of GM/c 2 being M = Mi + M 2 . 



MKN 421 

(M„M a )= (0.3,2.4)xl0 5 M e 
p, = 1 89x10 s 
e, = 0.5 




Fig. 4. The orbital evolution of the black hole binary at 
the center of Mkn 421 with initial parameters able to re- 
produce the observed lightcurve periodicity is shown. The 
orbital perihelion shift is evident. 



orbit will be tangent to a different Keplerian orbit. In the 
osculating orbit formalism a general two-body orbit is gen- 



erally described by six parameters: i the inclination of the 
orbit with respect to a reference plane, the line of the 
ascending node, to the angle between the line of node and 
the pericentric line, a the semi-major axis, e the orbital ec- 
ce ntricity and T the time of pericentric passage. According 
to lLincoln fc Willi 1)199(1) . these quantities are coupled by 
a set of first-order differential equations, see eq. (2.11a)- 
(2.11c) in the previously mentioned paper, which can be 
numerically solved, by appropriately choosing the initial 
conditions, obtaining a = a($), e = e($) and ui = w($), $ 
being the usual polar angle which in turn depends on time 
t through the relation r 2 <i> = (Afp) 1 / 2 . Here, the auxiliary 
parameter P is given by p = a(l — e 2 ). 

Applying this approach, we are able to describe the 
motion and the orbital evolution of the binary system. 
In Figs. |3| a) and 6) the evolution of the possible Mkn 
421 MBH binary system is shown. Notice that the or- 
bital parameters reported in the same panels correspond 
to a set of values of the Lorentz factor 7b, semi-major 
axis a and orbital eccentricity e reproducing the observed 
X-ray light curve periodicity. The evolution of both the 
MBH binary system eccentricity a) and separation b) are 
given as a function of the number of revolutions $/27r. 
In the sa me panels the expecte d coalescing time, calcu- 
lated as in lLincoln fc Willi l|l99Cl|) . is also reported both in 
geometrical units (2.6 x 10 8 ) and in years (2 x 10 5 yrs) 
for the considered MBH binary system. For comparison, 
the coalescing time scale of an eccentric binary cal culated 
using the Peters approximation (see eq. (5.14) in IPeters! 
l|l964|) ') gives ~ 1.3 x 10 5 yrs. It is not sur prising that the 
two coalescing time scales differ, since in iLincoln fc Wi3 
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Fig. 5. The two polarization states (h x and h+) of the gravitational waves emitted by the MBH binary system (with 
reduced mass /z), at distance R from Earth, for the selected model in Fig. [3] are shown. The waveforms are given for 
the first 10 orbital revolutions. The two parameters O and W determine the direction of observation. In principle, the 
study of the gravitational wave templates will allow one to extract information about the MBH orbital parameters. 



( 1990) different approximations and relativistic effects (as 
the perihelion shift) are taken into account. As evident, 
the difference between the time scales eval uated by usinf 

ilPetersI I 



these methods is not dramatic, so that the| 
formalism can be used, as a first approximation, to eval- 
uate the coalescing time for binaries that are (like ours) 
far enough from the plunge stage. A fur ther effect that 
can b e accurately described by using the ILincoln fc WiJ 
( 199oJ treatment is the orbital evolution of the black hole 
binary system. For example in Fig. 0] we show the or- 
bital evolution of the black hole binary in Mkn 421 with 
initial parameters chosen in order to reproduce the ob- 
served lightcurve periodicity. The orbital perihelion shift 
is clearly evident. 

Once the orbital motion of the binary system is known 
as a function of time, the osculating orbital parameter 
formalism allows us to determine the gravitational wave 
form, i.e. the evolution in time of the metr ic perturbation. 
In par ticular, using eqs. (4.1a)-(4.16) in ILincoln fc Wifl 
(1990), we can evaluate the polarization states h x and 
h+ of the gravitational wave emitted by the considered 
MBH binary system. In Fig. we plot, for the first 10 
orbital revolutions, the expected wave forms depending on 



the and \& parameters which determine the observati 
direction. 



4. Discussion and conclusions 

Usually, at the first stage of evolution, our binary systems 
are outside of the LISA frequency band since the typical 
frequency of the emitted gravitational radiation is much 
lower than 10~ 4 Hz. However, other exper iments, such as 
the A ST ROD gravitational wave detector l|Wu. Xu fc Ni I 
l200Ct In3|2002|) . will have a much higher sensitivity than 
LISA. In fact, even using the same laser power as for LISA, 
the ASTROD sensitivity would be shifted to a fre quenc y 
lower by a factor up t o 60 (30 on av erage), In3 <)2002l) . 
Moreover, as noted bv iRiidigerl l|2002fl . if the LISA ac- 
celerometer noise goal will be obtained, the ASTROD sen- 
sitivity at low frequencies will be about 30 times bette r 
than that of LISA, as indicated in Fig. 8 of [n] l|2002t) . 
4 Fig. El shows the possibility of detecting gravitational 
waves from such systems using LISA and ASTROD(l) and 
ASTROD (2) (the detailed descrip tion of these ASTROD 
facilities was given byE3 l)2002|) ). Moreover, if the ab- 
solute metrological accelerometer/intertial sensor can be 

4 See also Fig. 13 from the paper by IRiidige rl (2002). 
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Fig. 6. Assuming particular values of the free model pa- 
rameters (7b, e and a), we compare the expected grav- 
itational wave spectra with the LISA (L), ASTROD1 
(AD1) and ASTROD2 (ADF2) thresholds, respectively. 
Note that the values of the black holes reported for each 
Mkn objects strongly depend on the chosen parameters. 



developed, there is even the possibility of reaching the 
ASTROD(3) sensitivity curve with a shifting factor that 
could reach 10 3 - 10 4 lower than LISA ( see again Fig. 8 
from the paper by|NJ 1)200210 . Note that ways to decrease 
the instrum e ntal n oise up to (1 — 3) x 10~ 5 are discussed 
by iHughesI l|2002() (see also lLarson. Hiscock fe Hellingsl 
l|2002]) L 

In principle there is even a non-negligible chance to de- 
termine the inclination angle for a binary black hole sys- 
tem using gravitational wave observations if we will have 
a possibility to distinguish these templates for different 9 
and ^ angles. Of course, it could be only a hypothetical 
chance to extract this information from observations be- 
cause one should collect data for some years to reach the 
necessary sensitivity to detect the emitted gravitational 
waves. The dependence of gravitational wave templates 
on Q and \I/ angles could be important to construct opti- 
mal filters for gravitational wave detection. 
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